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Abstract 

Region-of-Interest (ROI) tomography aims at reconstructing a re¬ 
gion of interest C inside a body using only x-ray projections intersect¬ 
ing C with the goal to reduce overall radiation exposure when only a 
small specific region of the body needs to be examined. We consider 
x-ray acquisition from sources located on a smooth curve T in R 3 ver¬ 
ifying classical Tuy’s condition. In this situation, the non-trucated 
cone-beam transform Df of smooth densities / admits an explicit in¬ 
verse Z: however Z cannot directly reconstruct / from ROI-truncated 
projections. To deal with the ROI tomography problem, we introduce 
a novel reconstruction approach. For densities / in L°°(B ) where B 
is a bounded ball in R 3 , our method iterates an operator U combining 
ROI-truncated projections, inversion by the operator Z and appropri¬ 
ate regularization operators. Assuming only knowledge of projections 
corresponding to a spherical ROI C C B, given e > 0, we prove that if 
C is sufficiently large our iterative reconstruction algorithm converges 
uniformly to an e-accurate approximation of /, where the accuracy de¬ 
pends on the regularity of / quantified in the Sobolev norm W 5 (B). 

This result shows the existence of a critical ROI radius ensuring the 
convergence of the ROI reconstruction algorithm to e-accurate approx¬ 
imations of /. We numerically verified these theoretical results using 
simulated acquisition of ROI-truncated cone-beam projection data for 
multiple acquisition geometries. Numerical experiments indicate that 
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the critical ROI radius is fairly small with respect to the support re¬ 
gion B. 

Keywords', computed tomography, cone-beam transform, interior to¬ 
mography, region-of-interest tomography, ray transform. 


1 Introduction 

Computed Tomography (CT) is a non-invasive imaging technique, routinely 
used in medical diagnostics and interventional surgical procedures to visu¬ 
alize specific regions inside a body. CT involves patient exposure to x-ray 
radiation, with health risks of radiation-induced carcinogenesis which are es¬ 
sentially proportional to radiation exposure levels BE!- To reduce radiation 
exposure in CT, several strategies have been explored such as sparsifying the 
numbers of x-ray projections or truncating the projections so that only x-rays 
intersecting a small region-of-interest (ROI) are acquired. Reconstructing a 
density / from its projections is an ill-posed problem, meaning that small 
perturbations of the projections may lead to significant reconstruction errors. 
To address this problem, several approximate or regularized reconstruction 
formulas have been introduced over the years, such as the classical Filtered 
Back-Projection or the FDK algorithms [3j Ch.5]. However these methods 
are designed to work using non-truncated projection data. When projections 
are truncated, the reconstruction problem may become severely ill posed and 
non-uniquely solvable pE]. For instance, the so-called interior problem, where 
projection data are only known on a region strictly inside the support of the 
density /, has no unique solution in general [lj. As a result, naive numerical 
reconstruction algorithms such as direct application of a global reconstruc¬ 
tion formula, with the missing projection data set to zero, typically produce 
serious instability and unacceptable visual artifacts. 

The ROI reconstruction problem. The problem of ROI recontruction 
in CT has been studied in multiple papers and using a variety of methods 
(see, for example, the recent reviews m and the references therein). Recent 
remarkable results have shown that it is often possible to derive analytic ROI 
reconstruction formulas from truncated projections, provided the ROI is cho¬ 
sen with certain restrictions (cf. [7j E, ,3]). Such explicit ROI reconstruction 
formulas from truncated projections typically depend on the specific acquisi¬ 
tion modalities and impose restrictions on ROI geometry; for instance, some 
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prior partial knowledge of the density / within the ROI is required or the 
ROI cannot lie strictly inside the support of /. 

Iterative methods on the other hand provide a more flexible alternative for 
the reconstruction from truncated or incomplete projections as they can be 
applied to essentially any type of acquisition mode (cf. (TUI EH, EG])- Many 
such methods rely on total variation and other forms of regularization to 
ensure the convergence of the algorithm. For instance, the recent ROI re¬ 
construction approach by Klann et ah [13j relies on an appropriate wavelet 
based regularization. In this approach, the uniqueness of the interior prob¬ 
lem is guaranteed under the hypothesis that the density function is piecewise 
constant. However this result assumes the ideal case of a noiseless acquisi¬ 
tion and leaves the problem of stability in the presence of noise open. With 
respect to analytic formulas, iterative methods are usually computationally 
more intensive, especially for 3D data. However, advances in computational 
capabilities (e.g., [H]) and recent ideas from compressed sensing (e.g., [15] ) 
offer powerful tools to overcome this limitation. 

Our approach. In this paper, we consider the ROI reconstruction problem 
aiming at reconstructing an unknown density inside an ROI C using only 
the projections intersecting C. We fix a bounded ball B C M 3 and a smooth 
curve r C l 3 of ‘ray sources’ exterior to B. We assume that the objects 
illuminated by x-rays are strictly included in B and are characterized by 
their unknown density functions / G L°°(B). For any spherical region of 
interest C C B, we denote by 7 Zc the set of all half-lines (or ‘half-rays’) r 
emanating from arbitrary points of F and intersecting C. The C-truncated 
cone beam projection operator Dc maps any density function / e L°°(B) 
into a function Dcf defined, for each half-ray r e 7 Zc, by integrating / 
over r. 

The non truncated cone beam operator, corresponding to the case C = B, 
will be denoted D = Db- It is known that, if T and B verify classical 
geometric Tuy’s condition , Grangeat’s classical formula provides an inverse 
Z of the non truncated cone beam projection Z7, verifying ZDf = f for every 
C 2 -density / with compact support included in B. For several specific types 
of curves T, the non truncated operator D can classically be inverted by an 
operator Z defined on smooth densities by geometry-specific formulas. 

In the setting described above, we develop a method to construct approx¬ 
imate inverses Zq for the C-truncated operator Dc, defined by an iterative 
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RO/-reconstruction algorithm which converges whenever the unknown den¬ 
sity / is smooth enough and the volume of the difference set B\C is small 
enough. More precisely, we start with an explicit operator Z inverting the 
non-truncated cone beam projection D for smooth densities and we construct 
a regularization operator r such that U = Zt(D-Dc) becomes a contraction 
on L°°(B). Then, for any unknown density / in W 5 (B), we set / 0 = ZrDcf 
to define iteratively the approximating sequence of densities (fj) by 

/;+1 = fo + Uf„ (1) 

Our mathematical analysis proves that as j —>■ oo, the sequence fj converges 
to an e-accurate reconstruction / of /, at exponential speed in L°°(C), for 
any spherical ROI C C B having a radius larger than a critical radius p(e). 
That is, given an accuracy level e and a sufficiently large spherical region 
C C B, we generate an estimate / of / such that 

11/ - /IU°°(C) < e l|/||w5(S)- 

Our results also extend to the situation where sources are located over a 
whole sphere in M 3 containing B. 

Note that our AOJ-reconstruction approach can be applied whenever the 
non-truncated cone beam projection operator D can be inverted by an irn- 
plementable formula or a “blackbox algorithm” Z that is well defined on 
smooth densities. Unlike other methods proposed in the literature we do not 
need any explicit restriction on the ROI location or any prior knowledge of 
the density in the ROI as long as an inverse Z of the non-truncated cone 
beam projection D exists in the sense stated above. As indicated above, the 
existence of such e-accurate inversion of the truncated cone beam projection 
Dc in only guaranteed for C relatively close to B. 

The iterative scheme Q is formally similar to other iterative algorithms 
also proposed in the literature for ROI reconstruction such as the so-called 
Iteration Reconstruction-Reprojection (IRR) algorithm [TSj (T7, IS] , the Or¬ 
dered Subsets Convex algorithm [[T9] proposed to speed up CT reconstruction 
by reducing the number of projections and the iterative maximum likelihood 
(ML) algorithm proposed by Ziegler et ah |T8j However, existing applications 
of the IRR method and other iterative methods for ROI reconstruction found 
in the literature are mostly heuristic and provide no theoretical justification 
for convergence. In this paper, we provide a rigorous analysis of the inver¬ 
sion of the cone beam transform for sources located on a three-dimensional 
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curve satisfying classical Tuy’s condition. Using this theoretical framework, 
we prove that it is possible to define and compute an approximate inverse of 
the truncated cone beam transform. 

To validate our approach in the discrete setting, we have performed nu¬ 
merical experiments using four classical discrete x-ray acquisition geometries, 
with sources located on a sphere, a spiral, a circular curve and twin orthog¬ 
onal circles. For each setting, we have simulated ROI-truncated cone beam 
data acquisition using three different density functions in M 3 : a Shepp-Logan 
phantom, a mouse tissue density data sample, a human jaw density data sam¬ 
ple. We have performed extensive numerical tests using spherical ROIs with 
various centers and radii and found that the numerically computed ‘critical 
ROI radius’ is relatively small as compared to the size of the support of / 
and essentially insensitive to the ROI location. 

Paper outline The paper is organized as follows. In Section [2j we re¬ 
call the definitions of the ray and cone-beam transforms, and classical Tuy’s 
condition valid for acquisition settings with sources on smooth 3D curves. 
In Section [3j we examine known inverse operators Z implementing the re¬ 
construction of densities from non-truncated projection data and study the 
continuity properties of Z on adequate Sobolev spaces defined on the space 
of rays 7 Zb- In Section [4j we define a class of smoothing approximations of 
the identity in the image and projection domains, and we indicate how to im¬ 
plement these regularization operators by ‘small’ mollification. In Section [5j 
we describe our iterative ROI reconstruction algorithm from ROI-truncated 
data and we prove our main convergence results. In Section [6j we present 
numerical implementations of our iterative ROI reconstruction for discrete 
acquisition setups where sources are located on (1) a sphere, (2) a spiral, (3) 
a circular arm, (4) twin orthogonal circles, with simulated ROI-truncated 
x-ray data acquired from three densities in M 3 : a Shepp-Logan phantom, a 
mouse tissue density and a human jaw density. We analyze the accuracy of 
our ROI reconstruction approach and explore how the ROI radius impacts 
accuracy. Finally, we make some concluding remarks in Section [7} 

2 X-ray projections and Tuy’s condition 

We consider classical projection operators mapping density functions with 
domain in M 3 into linear projections defined on appropriates spaces of rays. 
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The most prominent examples of such projection operators are the ray trans¬ 
form and the cone-beam transform [3|. 

Recall that a ray r(u, 9) in R 3 is a line passing through the point «6l 3 
and parallel to the vector 9 G S 2 , where S 2 is the unit sphere of R 3 . That is 
r(u, 9) = {u +19 : t G R}. A half-ray r(a, 9) in R 3 is a half-line originating at 
the point a G R 3 and parallel to the vector 9 G S 2 . That is r(a, 9) = {a + t9 : 
t > 0}. 

2.1 The ray transform 

The ray transform maps a function / 6 L 1 (R 3 ) into its linear projections Xf 
obtained by integrating over rays at various locations and orientations, that 
is, 

/ OO 

f(u + t9 ) dt, 

-OO 

for mgR 3 and 9 G S 2 . Since Xf(u, 9) does not change if u is moved parallel 
to 9, it is sufficient to restrict u to the plane through the origin that is 
orthogonal to 9 in R 3 , henceforth denoted by T(9). Thus, Xf is a function 
on the tangent bundle of the sphere that we denote by 

T={(u,9) : 9 e S 2 ,u e T(9)}. 

Note that the pairs (u, 9) and ( u,—9 ) give the same ray f(u,9 ), so that 
the mapping (u, 9) —> f(u , 9) is a double covering of T which can thus be 
viewed as a 4-dimensional Riemannian quotient manifold. The associated 
Riemannian volume element on T is dudQ(9), where dQ{9) is the surface 
area on S 2 and du is the Lebesgue measure on the plane T(9). 

We will consider the action of the mapping X on functions with compact 
support inside a fixed open ball B C R 3 of radius p centred at the origin. 
We denote by 7 b the subset of T associated with the rays passing through 
B, that is 

T b = {(u,9 ) G T : r(u, 9) D B ^ 0}. 

Thus, 7 b is an open submanifold of T with compact closure in T and the 
natural Riemannian volume element at (w, 9) G 7b is given by dudQ{9 ). We 
denote as L p (Tb), 1 < P < oo, the standard L p function spaces associated to 
this Riemannian volume. 
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2.2 The cone-beam transform 

The cone-beam transform maps a function / e L 1 (M 3 ) into the function Df 
defined by 

poo 

Df(a,9) = / f(a + t6)dt , 

Jo 

for a£l 3 and 9 e S 2 . Here, we view a as the source of the half ray r(a, 6) 
with direction 6. Hence, Df is a function on the space of the half-rays 

n = {(a, 9) : 9 e S 2 , a E M 3 }. 

The space 1Z has the structure of a smooth 5-dimensional Riemannian man¬ 
ifold, with natural local coordinates defined by a 6 l 3 and standard spher¬ 
ical coordinates on S 2 . In particular 7 Z has a Riemannian volume element 
dp = dadQ(9), where dQ{9 ) is the surface area on S 2 and da is the Lebesgue 
measure on M 3 . 

We assume that all the unknown density functions / have compact sup¬ 
port inside a fix an open ball B C R 3 of radius p centered at the origin. In 
the more realistic tomographic setups considered below, sources are located 
on a smooth bounded curve Tel 3 supported outside the ball B. We denote 
by IZb the subset of all half-rays with sources on T and actually intersecting 
B that is 


n B = {(a, 9) e 1Z : r(a, 9) n B ± 0, 8 e S 2 , a G T}. (2) 

We call IZb the set of active rays. JZ B is a 3-dimensional manifold of class 
C°° with natural local coordinates defined by the arclength parametrization 
t of the curve T and the standard spherical coordinates on S' 2 . Thus, 1Z B 
is a submanifold of 1Z with Riemannian volume element at (a, 9) G TZ B 
given by dtd9), where dt is the Lebesgue measure on M. The total volume 
m(B ) = p{lZ B ) is clearly finite. 

Sources on a curve: Tuy’s condition. When the ray sources are located 
on a piecewise smooth curve T exterior to a bounded open ball B. classical 
Tuy ’s condition on T and B (see M) ensures that every smooth function / 
with compact support included in B can be recovered from its non-truncated 
cone-beam projection Df. 
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Definition 1. Let B be an open ball of finite radius centred at the origin. 
Let T ei 3 \5 be a C°° curve of length L parametrized by t —>• j(t) G M 3 , 
for 0 < t < L, with non zero velocities 7 '{t). 

T is said to verify strong Tuy ’s condition if there is a C 1 function A = 
A(x, 9) defined for ( x, 9) G B x S 2 and with values in [0, L\ such that, for all 
( x, 9) e B x S 2 , 

(9, 7 (A(x, 9)) = 0 and (9, 7 '(A(x, 9))) 0. (3) 

Note that, by the Implicit Function Theorem, the function A is of class C°°. 

The strong Tuy condition is satisfied, for instance, when T is a long 
enough circular helix ” containing” the ball B, or when T is the union of two 
concentric circles positioned on orthogonal planes in M 3 . 

Even though an helix does not necessarily satisfy strong Tuy’s condition 
(in general, there are planes that intersect a helix at one point, with tangen¬ 
tial intersection), a bounded circular helix is complete (in the sense of Tuy) 
as long as the support of the density is sufficiently small, and is surrounded 
by the helix. If this assumption holds, then tangential intersections between 
the planes and source curve are negligible as they occur on a set of Lebesgue 
measure zero [51]. 

As mentioned above, we will consider in Section [6] discrete applications of 
the cone-beam transform for different practical acquisition setups including 
the spherical case , where T is a sphere surrounding the target ball B , the 
spiral case, where T is a segment of circular helix, the C-arm case, where T 
is a circular arc and the twin orthogonal circles case, where T is composed 
of two concentric circles positioned on orthogonal planes in M 3 . In all these 
cases there is a formula to reconstruct a compactly supported smooth density 
function / from its non-truncated projections. 

3 Reconstruction from non-truncated projec¬ 
tions 

For the non-truncated projection operators considered above, which map 
density functions in M 3 into a full set of linear projections, it is possible in 
many classical cases to define a formal inverse operator. 

For the non truncated ray transform X, when / is in the space S of 
functions on M 3 having fast decreasing derivatives of all orders and when 


Xf(u, 6) is known for all (u, 9) G T, then there exists an inverse operator Z 
such that f(x) = Z.Xf(x) (cf. Sec. 2.2] or [22]). 

For the non truncated cone-bean transform D , if the source location Y is a 
piecewise C°° curve exterior to a ball B and verifying strong Tuy’s condition, 
then for all / in C 2 (R 3 ) with compact support included in B , there is an 
inverse operator Z such that Z.Df(x ) = f(x), and Z can be implemented 
by one of several variants of Grangeat’s formula [23] G]. For example, in 
spiral tomography , where Y is a segment of a circular helix, the inverse Z 
of the non truncated cone beam transform D can be computed either by a 
variant of Grangeat’s formula or alternatively by the Katsevitch’s formula 
[211 25] . In the setting of C-arm tomography , where Y is an arc of circle, an 
approximate inverse operator Z of the non truncated cone beam operator D 
can be computed using again a variant of Grangeat’s formula [201126] . 

We point out that all these exact formulas inverting the non truncated 
cone beam operator require smoothness conditions on / to reconstruct / 
as a function. As noted by Natterer [3], Tuy [20] and other authors, to 
define the most generic linear operator Z inverting the transform / —> Df , 
one should consider / and g = Df as distributions instead of functions. 
However, numerical reconstructions of / from discretized projection data Df 
usually smooth the non truncated data Df before reconstruction. Therefore 
classical proofs of exact inversion formulas for non truncated cone beam data 
tend to focus on smooth density functions /. Indeed, in spiral tomography, 
where T is an helix, the original proof of Katsevich’s inversion formula in [21] 
requires / G Cq°(B)] finite degree of smoothness can be achieved using more 
sophisticated arguments [27]. Similarly, when Y is a smooth curve, the proof 
of Grangeat’s inversion formula in [23] requires / G C 2 (B). Also in the more 
academic setting of the ray transform, where the full set of projections (for 
all (u, 6) G T) is known, the inversion formulas in [3] 22] require the Fourier 
transform of / to decrease rapidly at infinity. 

In the following, we will define exact inverses Z of the non-truncated 
projection operators D as explicit linear operators acting on Sobolev spaces 
of densities. This definition will be useful to derive important continuity 
properties of Z. 

We start by defining appropriate Banach spaces to handle the space of 
rays. 
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3.1 Banach spaces of smooth functions on manifolds 

Let At be a Riemannian manifold of class C°° with volume element d/i and 
finite volume ^(At). Here we consider only manifolds which are either com¬ 
pact or are the interior of a compact manifold with a C 1 -boundary. One 
can then find and fix a finite covering of At by open relatively compact sets 
Uj,j G J endowed with diffeomorphic local maps hj : Vj —* Uj, where the 
Vj of are bounded open balls in R 3 , and each hj is the restriction to Vj of 
a local map defined on an open Euclidean ball containing the closure of Vj. 
Explicit such finite coverings Uj are easily specified for the manifolds of rays 
IZb given by (|2]), and for T x S' 2 , where Y is either a piecewise C°° bounded 
curve in R 3 or a whole sphere in R 3 . 

Fix as above a finite covering Uj,j G J, of At and the local maps hj : 
Vj —>■ Uj. A function g on At is said to be uniformly bounded if and only 
is all the g o hj are bounded. For any r > 0, call C r ( At) the space of all 
functions on At having continuous and uniformly bounded differentials of all 
orders up to r. For 1 < p < +oo, we denote L v { At) the usual Banach spaces 
of functions g on At such that \g\ p is /i-integrable and /i is the Borel measure 
on At. 

For each r, the space C r { At) is included in the Sobolev space W r (M.) of 
functions g G L 2 (At) endowed with the Banach space norm 

\\g\\w r (M) = ° hj\\w r (Vj) 

j&j 

We have the following standard result (cf. [3, 3]). 

Proposition 1. Let Y C R 3 be either a smooth bounded curve or a full 
sphere exterior to a bounded open ball B. Let IZb be the manifold of half¬ 
rays with sources on Y and intersecting B. For each integer r > 0, the non- 
truncated cone-beam transform D is a bounded linear operator from W r (B) 
into W r (flZ B ), as well as from L°°(B ) into L°°{JZb), and maps C r (B ) into 
C r (lZ B ). The non-truncated ray transform X is also a bounded linear op¬ 
erator from W r (B ) into W t {Tb) as well as from L°°(B) into L°°{Tb), and 
maps C r (B) into C t (Tb )■ 

3.2 Inversion of the non-truncated ray transform 

An analytic inversion formula for the ray transform can be derived from the 
classical Fourier slice theorem. In this section, we derive explicit continuity 
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properties for this inversion formula. 

For any 9 in S 2 , let T{9) be the plane orthogonal to 9 in M 3 and containing 
the origin. As seen in Sec. 2.1, any ray r(u,9) is non-ambiguously indexed 
by 9 G S 2 and u G T{9) and the tangent bundle T of the unit sphere in M 3 
is a 4-dimensional manifold with volume element dud9). 

Fix a an open ball B of radius p in M 3 and let Tb be the manifold of all 
rays intersecting B. For any function g(u,9 ) on Tb, let go be the function 
defined on T(9 ) by g$(u) = g{u,9). For v G Tg, the 2-dimensional Fourier 
transform of go on the tangent plane T{9) is given by 


T e g e (v)= [ e- i ^g(u,9)du. (4) 

Jt(6) 

whenever the integral is well-defined. Using standard inequalities, a direct 
computation shows that for g G IF 4 (7b) then, for any 9 G S 2 and u G T(9), 

\T e ge(v)\ < c(l + |u| 4 ) _1 ||^|| H /4 {Ts) , (5) 

where the constant c depends only on p and not on g. 

For / G L 2 (B), the usual 3-dimensional Fourier transform of / will be 
denoted by 


f(z)—Tf(z) — f e 1< ' Z,X ' > f(x)dx, for z G M 3 . 

J B 

By the Fourier slice theorem (cf. [3| Sec. 2.2]), for any 9 G S 2 and z G M 3 
such that (z, 9)= 0, the ray transform g = Xf of / verifies 

f( z ) = F9e(z), (6) 

provided the two Fourier transforms involved in the formula are well defined. 
As shown in [3J 22J when f(z) tends to zero at infinity faster than any poly¬ 
nomial in z, then equation ([6]) can be used to derive inversion formulas to 
reconstruct / from its non-truncated projections g. 

For the non-truncated ray transform X, we now specify a bounded linear 
inverse defined on W a (Tb)- A function g defined on Tb can be extended to 
T by setting g — 0 on T \ Tb 

Proposition 2. Fix a ball B and define Tb as above. Fix any Borel mea¬ 
surable function z —» 9(z) from M 3 to S 2 such that ( z,9(z )) = 0 for almost 
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all z G M 3 . For any g G W 4 (Tb) and, all x G M 3 , the following integral is 
necessarily finite: 

Jg{x) = (2vr)- 3 [ e^ Tg e{z) (z) dz. (7) 

Jr 3 

The restriction Zg = 1 B Jg of Jg to the ball B defines then a bounded linear 
operator Z from W 4 (Tb) into L°°(B) and from W 4 {Tb) into L 2 (B). More¬ 
over, for any f G W 4 (B), the non-truncated ray transform g = Xf verifies 
the identity f = ZX f. 

Proof. For any g G IF 4 (7b), the inequality ([5]) holds for all 6 G S 2 and v in 
the tangent plane T{6). Hence for all iGl 3 , the integral Jg(x), defined by 
equation (JtJ) , is bounded by 

\Jg(x)\ < c \\g\\ W i(j B ) [ (l + H 4 ) -1 ^ 

Jr 3 

< c\\g\\w±{T B )i 

where the constant c (changing from line to line) depends only on the radius 
of B. 

It follows that, for all g G W a {Tb ), there is a new constant c such that 
the function Zg = 1 b Jg verifies 

\\Zg\\L™{B) < c\\g\\wi(T B )- ( 8 ) 

Then g —)■ Zg is a bounded linear operator from W a (Tb ) into L°°(B) and 
hence also into L 2 (B). Moreover, for any / G W A (B), the function g = Xf is 
in W 4 (7b). Therefore the Fourier slice formula Q combined with ([7]) show 
that 

/ = 1 B jg = Zg = ZXf. 

This achieves the proof. □ 

3.3 Inversion of the non-truncated cone-beam trans¬ 
form 

We now construct an operator inverting the non truncated cone-beam trans¬ 
form D when the projections belong to a Sobolev space of rays. 

Fix a ball B and define Tb as above. Let T C I 3 be a C°° curve with 
support exterior to the open ball B and parametrized by 7 : [0, L] -> K 3 . 
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Assume strong Tuy’s condition is verified, i.e., there is a C l function A = 
\{x,6) :BxS 2 ^ [0, L\ verifying (§). 

Consider any function g G C 2 (1Zb) with compact support inside 7 Zb- We 
extend g to T x S 2 by setting g = 0 on V x S 2 \ TZb and then extend g to a 
C 2 function G defined on T x M 3 by 

G(s,») = || 9 ||- I 9 (s,/ iT ), forallser, B eM 3 \{0}. (9) 

IMI 

Now, for all t G [0, L\, y G M 3 \ {0} we set 

K (t, y ) = ^ d yG(l(t),y) (10) 

and, hence, for all x G B , we dehne the function Zg by 


Zg(x ) 


1 

87r 2 


>a&T(9) 



(9,K(X(x,9),a)) 


ddda , 


( 11 ) 


where, as above, T(0) C S' 2 is the set of all a G S 2 such that (a, 9) = 0. 

We have the following result. 

Propositions. The Grangeatformula 0 defines a linear operator g —> Zg 
from C 2 (1 Zb ) into L°°(B). Moreover there is a constant c depending only on 
T and the radius of B such that, for all g G C 2 (1Zb) with finite Sobolev norm 
\\9\\wHKb)> we have 

\\Zg\\L<x>{B) < c|fy||w 4 (^s)- (12) 

In particular Z can be extended to a bounded linear operator from W 4 (TZb) 
into L°°(B ) such that whenever g = Df is the non truncated cone-beam 
transform of f G W a {JZb), one has the identity f = Zg = ZDf. 


Proof. When g = Df with / G C 2 (B), the assertion Zg = f is proved with 
different notations in [3j Sec. 5.5.2] using a variant of the Grangeat’s inversion 
formula due to Zeng, Clack and Gullberg 


y 


For a generic g in C 2 (JZb), the vector valued function K, given by (10), is 
continuous by construction and hence remains bounded in M 3 for t G [0, L], 
eS 2 . 

In the following, c, C\, C 2 ,... denote positive constants which depend only 
on the radius of B and T but not on g. 

In equation ( [IT] ), the denominator den(x, 9) = (9, 9))) is continuous 

for x G B, 9 G S c and is never zero due to Tuy’s conditions, so that |den| > 
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c > 0 for some constant c. Then equation (11) readily provides a constant 
ci such that 

\\Zg\\ L °°(B)<ci sup \\K (t,a)|| R 3 . (13) 

t£[0,L\,a£S 2 

Set h(t,9) = g{lit)i9)- Equations ([9]) and (10) show that there is a constant 
c 2 such that for all g in C 2 (1Zb), 


\\K (t,a )|| R 3 < c 2 sup \\^-d e h(t,9)\\ R 3, 
te[o,A],6eS 2 “t 


(14) 


for all t G [0,L],a G S 2 . By definition of W a {JZb ), there is a constant C 3 
such that for any function g in W 4 (JZb) the function h verifies 


II^11 W 4 ([0, J 4]xS 2 ) < C3 |M|lV 4 (ft B )- 


(15) 


The Sobolev imbedding theorem in dimension 3 holds on the Riemannian 
manifold 7 Zb, relating the norm of h in C 2 (1Zb) with its norm in W a (1Zb)-> 
as explained in Appendix [7j It thus provides a constant c 4 such that, for any 
function h G W A ([0,L] x S' 2 ), all partial differentials j^dgh of order < 2 of h 
are bounded and continuous on [0, L] x S 2 and verify 


sup 

(t,6)£[0,L]xS 2 


-mi, 6) 


< C 4 ||^||vK 4 ([0,L]x5 2 )- 


(16) 


Combining the inequalities (13) (14) (15) (16), we get 


\\Zg\\L°°(B) < CiC 2 C 3 C411 (711 iy4 {JZb) 


which achieves the proof. 


□ 


Remark: Katsevich’s inversion formula. As mentioned above, in spi¬ 
ral tomography, the source curve T C I 3 is a circular helix, which can be 
parametrized as 7 (t) = (acos(t), asin(f), kt ) for some fixed positive a, k. Let 
/ be a compactly supported density function / G C™(B), with p < a (so that 
the helix is surrounding the support of /). Katsevich proved (cf. [29]) that, 
in this setting, / can be reconstructed from its non-truncated cone-beam 
projections Df by a formula which can be written as 

f = \(U + V).Df, 
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where U and V are explicit operators involving divergent integrals. This 
divergence is carefully analyzed by adequate approximations in [29], but this 
inversion formula remains rather unwieldy and the continuity properties of 
(U + V) are not easy to evaluate directly. 

As mentioned above, even though an helix does not satisfy strong Tuy’s 
condition in general, this condition holds for a bounded circular helix as long 
as the support of the density is sufficiently small. Hence, we can define an 
explicit inverse operator Z of the non-truncated operator D according to 
Proposition [3j With respect to Katsevich’s inversion formula, our approach 
has the advantage of reconstructing / from projection data Df G W a {JZb) 
for all / G bP 4 (i?) and to provide a precise Sobolev continuity property for 
the inverse operator Z. 

4 Regularization in the space of rays 

We now define and construct a class of regularization operators in the space 
L 2 {1Zb)- We start by defining a notion of approximations of the identity on 
the manifolds considered in Section [3] 

Definition 2. Let Ai be a C°° Riemannian manifold with volume element 
dp and finite volume. As above, assume that Ai is either compact or is the 
interior of a compact manifold with a C 1 -boundary. For any integer r > 1, 
we call C r approximation of the identity in L 2 (Ai,p) any sequence of linear 
operators tn ■ L 2 (Ai) —> C r (Ai) verifying the following conditions. 

(i) There is a constant c such that for all g G L 2 (Ai) and all integers N 

\\TNg\\wr(M) < C-N r \\g\\ L 2(M)- (17) 

(ii) For any g G L 2 (M ) 

lim || g - T N g\\ L 2 {M) = 0, for each g G L 2 (M). 

V->oo 

(Hi) For each integer 2 < p < (r + 1) there is a constant c such that for all 
g G W p (M) 

c 

\\g — T Ng\\wp- i (M) < jj\\g\\wp(M)- 

(iv) Whenever g has compact support then r^g also has compact support. 


15 


Note that N —)• tjv will remain a C r approximation of the identity in 
L 2 ( M.,v) for any positive Borel measure v on A4 such that both densities 
^ and are bounded. 

afi av 

We next show how to construct C r approximations of the identity in 


4.1 Approximation of the identity by small convolu¬ 
tions 


When the manifold At is a bounded open Euclidean ball in M fc , one can 
generate a C r approximation of the identity as follows. Select any fixed C r 
function w > 0 on with compact support and Lebesgue integral equal to 
1 and, for / G L 2 (M fc ), define the “small” convolutions by 


cn! = / * %, 


where wn(x) = N k w(Nx). Standard results on convolutions show that the 
sequence ( a n) verifies the properties (i),(ii) and (iv) of Definition [2j The 
proof of property (iii) is the following. 

For z G and 1 < p < oo, we define 4> p (z) = (1 + \z \ 2 ) p ^ 2 . Denoting by 
Wn the Fourier transform of Wn, we have that Wn(z) = w(z/N), where w 
is the Fourier transform of w. Hence, for all z G M fc , 

\ Wn (z)-1\ = \w(z/N)-1\<^\z\, (18) 


where c 
that 


w 


\L°° ■ 


From inequality (18), since (ctjv/) = Wjv/, it follows 


(f)p-i(z) | (a N f)(z) - f(z )| = (t> P -i{z) | W N (z)f(z) - f(z )| 

< ~j^\ z \ \K Z )\ 

< jjM z )\f( z )\i 

for all z G M fc . From the last inequality and by the definition of Sobolev 
norms, we then get: 

||<Tv/ — fllwp- 1 ^) — W&P-I ((T/v/) — /)||,L 2 (R fe ) 

c * 

- ^H0p/IU 2 (K fc ) 

— jy II /II W p (R fe )• 
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This proves property (iii). □ 

Small convolutions are easily and explicitly extended to the manifolds of 
active rays 7 Zb, as we now show by patching together local small convolutions 
through appropriate local maps. 

Proposition 4. Let M. be a C°° Riemannian manifold, of finite volume. As 
above, assume that Ai is either compact or is the interior of a compact man¬ 
ifold with C 1 -boundary. Then, for any integer r, one can construct explicitly 
a C r approximation of the identity N —>■ on L 2 (Ad). 

Proof. As observed above, on Ai we can select an open finite covering Uj,j E 
J, and local maps hj : Vj -E Uj to construct a finite partition of unity by C°° 
functions Uj with compact supports included in Uj and verifying 0 < Uj < 1 
and Y2jej u j = 1- On each Euclidean ball Vj, select a C r approximation 
of the identity N —y (j ) in L 2 (Vj), for instance by small convolutions as 
indicated above. 

For any g E L 2 (A4), let gj = guj and define the operators g —> Tug by 


tn9 = Ll G i° h i '' 


(19) 


where Gj = cr A r(j)(gj o hj) Since each map hj can be smoothly extended to a 
neighborhood of Vj, each hj has bounded derivatives of any order. Hence the 
mapping gj -E gj o hj is a bounded linear operator from L 2 {Uj, /i) to L 2 (Vj), 
and from W r (Uj) to W r (Vj). Similar boundedness properties hold for the 
linear operators mapping Gj —$■ Gj o hj 1 and g —» guj. The explicit formula 
(19) and the fact that the (j ) are C r approximations of the identity in 


L 2 (Vj) then implies directly that the sequence (rjv) satisfies the properties 
(i)-(iv) of Definition [2] Hence (r?v) is a C r approximation of the identity in 
L 2 (M,p). D 

Proposition [4] applies in particular to manifolds of active rays 7 Zb associ¬ 
ated to an open ball B and a smooth set of sources T exterior to B. 


5 Reconstruction from ROI-truncated projec¬ 
tions 

Fix a bounded open ball B C M 3 and a smooth set of sources T verifying 
strong Tuy’s condition. Let C be a spherical ROI strictly included in B. As 
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Figure 1: ROI-truncated cone-beam acquisition: projections are restricted 
to half-rays intersecting the ROI, which is a ball C included in the target 

ball B. 

illustrated in Figure [lj the C-truncated cone-beam transform Dcf of / is the 
restriction of Df to the manifold IZc C 71b of half-rays which intersect C. 
Denoting the indicator function of a set G by 1 g, one can then write 

Dc = 1 n c D, Yc = l-jig-KcD an d D = Dc + Yo- (20) 

by Proposition [lj Dc is a bounded linear operator from L°°(C ) into L°°(1Zc)- 
As observed above, the non-truncated operator D can be inverted by a 
linear operator Z such that / = ZDf for all densities / G W 4 (B) having 
compact support. However the C-truncated operator Dc cannot be directly 
inverted by applying Z to Dcf, even within the region C. So we now formally 
define ‘approximate’ inverses Zc for Dc ■ 
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5.1 Approximate inverses of ROI-truncated cone beam 
transforms 

Definition 3. Fix a bounded open ball Bel 3 and a smooth set of sources 
T verifying strong Tuy’s condition. Let C be a spherical ROI strictly included 
in B. For any e > 0, we say that the ROI-truncated cone-beam transform 
Dc admits an e-accurate inverse Zc if Zq is a bounded linear operator from 
L°°(TZc ) to L°°(B ) verifying 

||(1 — Z c D c )f\\L™(B) < e||/||vr 5 (B) 


for all f G W\B). 

Note that the e-accurate inverse Zq of Dc is in general not unique. Indeed 
recall that there are constants c and C\ such that for any / G W 5 (B ) one has 
||/IU”(B) < c||/||iv 5( B ) and hence 

\\Dcf\\L°°(n B ) < c i ll/IU°°(B) < cic\\(\\f)w 5 (B)- 

Fix now any bounded linear mapping K from L°°(TZb) into L°°(B) with 
operator norm \\K\\ < For any / G W 5 (B) we then have 

\\kD c f\\L°°(B) < e||/||vn 5 (B)- 

For any e-accurate inverse of Dc the operator Zc + K will then be a 
2e-accurate inverse of D c . 

5.2 Reconstruction from ROI-truncated projection data 

Let B C l 3 be a spherical region and F inside B a smooth curve verifying 
strong Tuy’s condition. As seen above, the non truncated cone-beam trans¬ 
form D with sources on T has then an inverse Z : W 4 (Hb) — > L°°(B) such 
that ZDf = f for all / in W 4 (B). The theorem below shows how to obtain 
an e accurate inverse of the ROI truncated operator De- 

Theorem 1. (ROI-reconstruction algorithm). With the notation above, let 
N —>■ ttv be any C A -approximation of the identity in L 2 (JZ B ) as in Definition 
[£| For any sphere C C B define the operators 

Un = Ztn(D — D c )■ 
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Given any e > 0, one can then find N = N(e) and 77 (e) > 0 such that 
the operator U = LOv(e) becomes a contraction of L°°(B) provided rad(B) — 
rad{C) < 77 (e). 

With N fixed as above, for any g e define the functions fj j = 

1,2,... by the recurrence 

f j+ i = fo + Uf j (21) 

with /o = Zrjyg. Then as j —» 00 the sequence (fj) converge at exponential 
speed in L°°(B ) to a limit Zcg■ This defines a bounded linear operator Zc 
from L°°(IZ B ) into L°°(B). Moreover for any unknown density f in W b (B), 
one can use the ROI-truncated data g = Dcf to compute an approximation 
of f given by f = Z c D c f and verifying 

\\f - f\\L°°(B) < e\\f\\ws(B) ( 22 ) 

Hence Zc is an e accurate inverse of the ROI truncated operator Dc- 


Remark. According to Theorem [lj the truncation region C is strictly in¬ 
cluded in B and must be large enough for the theorem to hold. The proof 
of theorem (presented below) provides an explicit lower bound rad(R) — 
rad(C)? 7 (e) implying the existence of an e-accurate inverse for the C-truncated 
cone beam projector Dq. Our theoretical estimate of 77 (e) is clearly too pes¬ 
simistic. Our numerical tests (see Sec. |b| indicate for instance that, for 
e = 0.10 and for all spheres C C B with same center as B and radius 
rad (C) > rad(B)/ 2, our i?OJ-reconstruction algorithm (21) converges to an 
e-accurate inverse Zc of Dc- This is a favorable situation for radiation expo¬ 
sure reduction by ROI truncated data acquisition combined with our ROI 
reconstruction algorithm. 

As seen above, e-accurate inverses are generally not unique. The following 
corollary outlines alternative constructions of Zc- 

Corollary 1. The notations and hypotheses are the same as Theorem [IJ Fix 
two C A -approximations of the identity: (r n ) in L 2 (IZ B ) and (a n ) in L 2 (B). 
For any n e N and any spherical region C C B, define the operator U n = 
a n ZT n (D — Dc ). Then, given e > 0, one can find N(e ) and 77 (e) > 0 such 
that provided rad(B) — rad(C ) < 77 (e), the operator Un^ is a contraction of 
L°°(B). Then for each g e L°°(B), the sequence (fj), given by the recurrence 
(21) with U = UN(e) and /o = er nZt^, converges in L°°(B) to a limit Zcg, 
where Zc is an e-accurate inverse of the ROI-truncated transform Dc- 
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5.2.1 Proof of the Theorem [Q 

Before presenting the proofs, we need the following two lemmata. For the 
remaining of this section, let T, B , C, D, Dc be given as above. 

Lemma 1. There is a constant c determined by T and B only such that, for 
any sphere C C B and any f e L°°(B), the linear operator Yq — D — Dc 
verifies 

\\Ycfh^ B , < c(rad(B) - ra<i(C)) 1/2 ||/||i~( B ), (23) 

where rad(C) is the radius of C. 

Proof: Let s be any source position on the curve T. Denote by z(C) the 
center of C. Call H(s, C ) C S 2 the set of all 6 e S 2 such that the half-ray 
r(s, 6) intersects C. The set of all these half-rays is a cone of revolution with 
vertex s, axis [s, z(C)], and half-aperture angle 0 < a(s,C) < 7t/ 2. The area 
of the spherical cap H(s, C) is hence classically given by 

area (H(s, C )) = 27t( 1 — cos(a(s, C )). (24) 

Elementary geometry yields 

cos(a(s, C )) = fc(rad(C), |s — z(C) |), (25) 

where k(u,v) = (tt2+ * a)1/2 ■ 

Since s e T and C C B, the numbers u = rad(C) and v = |s — z(C)\ 
remain respectively in the bounded intervals [0, rad(5)) and [m,M], where 

0 < m — —rad (B) + min Is — z(B) I and M = rad(5) + max Is — z(B) |. 

ser 1 ser 1 

The function (u, v ) —» k(u, v ) is C°° on the rectangle J = [0, rad (B)] x [m, M]. 
Hence there is a Lipschitz constant c such that, for ah (ui,Ui) and (u- 2 , V 2 ) in 
J, one has 


\k(u h vi) - k(u 2 ,v 2 )| < c(|wi - u- 2 1 + \vi - v 2 \). 


Due to equation (25), this implies 

| cos(a(s, C)) — cos(a(s, B))\ < c(rad(H) — rad(C)) + |(|s — z(C)\ — |s — z(H)|)|. 

Since |(|s — z(C)\ — |s — z(B)\)\ < \z(C) — z(B)\ < rad (B) —rad(C), it follows 
that 

| cos(a(s, C)) — cos(a(s, B))\ < 2c (rad(H) — rad(C')). (26) 
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The total volume of lie is 


vol (lie) = I area (H(s,C))ds, 


where dsdO is the volume element in the manifold of active rays lie■ Hence, 


due to (24), we have 


vo1(7£b — lie) = / (area (H(s,B)^ area(H(s,C)) ds 


= 2n J (cos(a(s, C)) — cos(a(s, B)) ds. 


From the last equation, using (26), we obtain 


vol(7 Zb — He ) < 47tLc (rad(H) — rad(C')), 


(27) 


where L is the length of F. The L 2 norm of the indicator function 
then verifies 


WItZb-ticWI^tZb) = v °1(^b - K c ) < 47rLc(rad(H) - rad(C')). 

Due to Proposition [TJ there is a constant c 0 such that, for any / G L°°(B), 
one has ||Z7/|| Z/ oo( 7 ^ fl ) < c 0 ||/||i,oo(b)- Hence, since C C B we have 


\\Y C f\\L\Tl B ) ~ 
< 

< 


\\^n B -n c D 

11 ^-n B -n c 11 l 2 {n B ) 11 D f 11 L oo (B) 
ci(rad(H) - rad(C')) 1/2 ||/|| L ^( B ), 


where Ci = Co^tcLc ) 1 ^ 2 . This achieves the proof when T is a curve of length 
L and does not intersect the closure of B. □ 

As seen above, the non-truncated cone beam transform D can be in¬ 
verted through a bounded linear operator Z : W^ilis) —» L°°(B ) given by 
the Grangeat’s formula (11). The lemma below shows how to construct a 
contraction on L°°(B) from Z and D — Dc- 


Lemma 2. Let (tn) be a C 4 approximation of the identity in L 2 (Hb) as 
in Definition [2| Then there is a constant k such that for any C C B with 
radius verifying rad(B) — rad(C) < k/N 8 , the linear operator Ujy = Zt n Y c = 
Ztn(D — Dc) is a contraction from L°°(B) into L°°(B), with operator norm 
II U n\ \ L°° (B) < 0.9. 
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Proof. By inequality (23), there is a constant c such that, for all C C B and 
all / G the operator Yc = D — Dc verifies 

\\Ycf\\i?yi B ) < c(rad(S) -rad(C)) 1/2 ||/|| L oc (i?) . 


Applying inequality (17) to g — Ycf with r = 4, we obtain a new constant 
Ci such that, for all C C B, all / G L°°(B) and all N 

IItw^c/H w 4 (7Js) — c^WYcfWv^ 

< CicN A (rad(B) - rad(C')) 1/2 \\f\\ L °°(B)- 


By applying inequality (12) to the function r^Ycf , we then obtain a new 
constant c 2 such that, for all C G B, all / G L°°(B) and all N , 


||£^v/||l°°(b) — 
< 

< 


||^jV^C , /|U°°(B) 

C 2 || vK 4 (7e s ) 

c 2 cidV 4 (rad(5) - rad(C')) 1/2 1|/1| Z/00(jB) . 


Set c 3 = c 2 cic. Then, provided rad(S) — rad(C') < c 3 /N 8 , the linear oper¬ 
ator U N is a, contraction of L°°(B ) with operator norm ||Thv||L°°(B) < 0.9. □ 


We can now prove Theorem |TJ 

Select and fix a C 4 approximation of the identity N —» tjv) in L 2 (1Zb), 
so that the bounded operators '■ L 2 {JZb) —> W 4 (TZb) verify definition [ 2 ] 
We will use the following shorthand notations for the norms of various linear 
operators T : 


\T\w5W5 

is the norm of T : W 5 (B ) ->• W\lZ B ) 

\T\w5W4 

is the norm of T : W 5 (1Zb, 

) -)• W 4 {1Zb 

\T\wALoo 

is the norm of T : W 4 (1Zb, 

) —>■ L°°(5) 

\T L2W4 

is the norm of T : L 2 (1Z B ) 

->• ^ 4 (77 b ) 


Due to Propositions [T] and [3j Definition [ 2 ] and equation ( fll| ), there is a 
constant c > 0 such that for all integers N 

\D\w5W5 < C , |^| W4Loo < C, 1 I — Tn\w5W4: < c/N, \tn\l2W a < cN 4 . (28) 

Given an e > 0, £x an integer N — N(e) by 

N = N{e) = 10c 3 /e. (29) 
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Since N = N(e) is now fixed, we will write r = r^. By Lemma [2j there is 
a constant k such that for any spherical region C C B verifying rad(R) — 
rad(C) < k/N(e) 8 , the operator Un — U — ZtYc is a contraction from 
L°°(B ) into L°°(B), with operator norm ||t/||z,“(B) < 0.9. 

We now fix 77 (e) = k/N(e) 8 and assume that the ROI radius verifies 
rad(.B) — rad(C) < 77 (e), which forces U to be a contraction with norm 
inferior to 0.9. Given any g in L°°(1Zb), define a sequence (fj) C L°°(B ) by 
the algorithm (21). By construction one has then 


\\fj+l ~ fj ||l°°(B) < ll^llll/? — /j —11| < 0.9||/j — fj-i\\L°°(B)- 

Hence, as j — > 00 , the sequence (fj) converges at exponential speed in L°°(B ) 


to a limit Zcg G L°°(B). By (21), we must have 


Zc9 = fo + U Z c f = Zrg + U Z c g 


(30) 


This defines a linear operator Zc '■ L°°(1 Zb) —> L°°(B). We now show that 
Zc has bounded operator norm. Since ||G||l°°(_b) < 0.9, the operator (I — 
U ) : L°°(B) —> L°°(B) has a bounded inverse given by the converging series 
£~ 0 U j , which yields 


II {i-u) 


-11 


|l°°(s) < 


1-0.9 


= 10 . 


Equation (30) yields that, for all g G L°°(1Zb), 

Z c g = (/ - Uy'Zrg. 


(31) 


Hence, due to the bounds (28), 


\\Zcg\\L°°(B) < || (I-U) 1 \\l°°{b)\Z\ WALoo | T IL2 W4 11 g \ \ L 2 (U B ) 

< 10c 2 iV( e ) 4 ||(7|| i2(7 , B) 

< lQc 2 N(eYm(B) 1 / 2 \\g\\ L c 0{B) (32) 


where m(B) is the finite Riemannian volume of TZb- So the operator norm 
of Zc is bounded. 

For / G W 5 (B) we have / = ZDf and the ROI truncated transform 
h = D c f of / belongs to W 5 (IZb) C L°°(IZb)- We then can write 


/ — ZrDf = Z (I — r)Df. 
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Combining this observation with the bounds given by (28) for the norms of 
Z, (/ — r) and D, we obtain, for all / £ W 5 (B), 


\\f - ZrDf\\ L ~ {B) = ||Z(/ — r)Df\\ L oo( B) 

„3 


< 


N(e) 


\\f\\w 5 (B)- 


Since D = D c + Yc, for any / £ W 5 (B) we have 

/ - ZrDf = f- ZrYcf - ZrD c f = (J - U)f - Zrh. 


(33) 


(34) 


From (J33J) and ( |34| we then get that for all / £ W 5 (B ) 

W ~ U)f ~ Zrh\\L°°( B ) < J~jrj—j\\f\\w 5 (B) 

and, hence, since || (/ — U)^ 1 ||l°°(_b) < 10 , we conclude that 
||/-(/-C)- 1 Zrh|| i oo (B) < ||(/ — C) - 1 ||z / oo ( £ ) ||(J — U)f — ZtIi\\l°°( B ) 


< 10 


N(e, 


Mf\\w 5 (B)- 


(35) 


For any / £ W 5 (B), the expression of Z c h = Z c D c f given by equation (31) 
then implies 


/ - Z c D c f = / — Z c h = / — (/ — 


Hence equation (35) yields 


11/ - ZcDcf \\l°°(B) < 10-^y ll/ll |w 5 (B)- 


Our choice of TV in (29) forces 10 
regions C C B verifying rad(H) — 


= e. Thus, for all / £ W 5 (B ) and all 
rad(C) < 77 (e), we get 


II (I ~ ZcD c )f\\L°°(B) < e||/||iv 5 (B)- (36) 

By definition |3j Zc is thus an e-accurate inverse of Dc- This completes the 
proof of Theorem [lj □ 


25 












Proof of Corollary [l] This proof is similar to the argument used in the 
proof above and will just be sketched. Select two C 4 approximations of the 
identity (crjv) in L 2 (B) and (tjv) in L 2 (JZb) and set Un = VnZtnYc. By 
applying Lemma [2 to Un we have that, given any e > 0, there is N — N(e) 
large enough and 77 (e) small enough to ensure that, for all regions C verifying 
rad(A) — rad(C) < 77 (e), the operator U = Un satisfies ||t/||z,°°(,B) < 0.9. The 
iterative sequence f 3+ \ = f 0 + Uf 3 initialized by the new f 0 = aNZr N g will 
then converge to a limit Zcg in L°°(A). As in the argument above, this 
defines an e-accurate inverse Zc of Dc■ □ 


5.3 Extensions 


Theorem [l] and Corollary [T| can be easily extended to the case where the 
set of sources T is a sphere strictly including B rather than a curve, using 


arguments nearly identical to those used above. The constant c in (23) then 


depends on the surface of the sphere T instead of the length of the curve T. 

Our RO I-reconstruction algorithm to approximately invert the AO/-truncated 
cone beam transforms Dq can also be extended using nearly identical argu¬ 
ments to generate approximate inverses of the AOJ-truncated ray trans¬ 
forms Xc. 

For discrete data Dcf generated by a numerical AO/-truncated cone 
beam transform of the unknown density /, our ROI reconstruction algo¬ 


rithm (21) can formally be applied, starting with any discretized reconstruc¬ 


tion algorithm Z specific to the acquisition geometry at hand and known to 
perform well on non-truncated data. From a formal point of view, Z can even 
be a ‘black-box’ inversion software dedicated to inversion of non truncated 
data. Of course the proofs Theorem [T| do not a priori cover such cases from 
a theoretical point of view. Yet, our extensive numerical tests in Section [6] 
indicate that our ROI-reconstruction approach does perform well in many 
practical CT setups. 


6 Numerical experiments 

In this section, we present extensive numerical experiments to evaluate the 
performance of our ROI reconstruction algorithm in multiple discrete settings 
with sources located on a smooth curve or a sphere. Since we had no access 
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to ROI-truncated data acquired with an actual CT device enabling cone- 
beam ROI truncation, our numerical experiments simulated ROI-truncated 
cone-beam acquisition. We used four classical acquisition geometries and 
multiple spherical ROIs with different locations and sizes applied to several 
3D density data. The goal of these numerical experiments was to numerically 
quantify accuracy of our ROI reconstruction algorithm within the ROI and 
to investigate how the ROI radius impacts this measure of accuracy. 

For a given cone-beam acquisition setup with target ball B C M 3 , a 
point z in B and a number v > 0. our Theorem [l] and Corollary [l] imply the 
existence of a critical radius p such that, for any spherical region C C B with 
center z and radius rad(C) > p, and for any density / with ||/||z,°°(s) < v, 
our ROI reconstruction algorithm from truncated data will converge within 
C to a good approximation of the unknown /. Our numerical experiments 
provide practical evaluations of this critical radius p. 

6.1 Simulations of ROI-truncated acquisition 

We have simulated ROI-truncated cone-beam acquisition for three discretized 
densities /: a 3D Shepp-Logan phantom; a 3D scan of mouse tissue; a 3D 
scan of a human jaw. For each discretized density /, given by a 3D image 
of size 256 3 voxels, we have first computed discrete non trucated cone-beam 
projections Df by simulating discrete acquisition and used these data to gen¬ 
erate ROI-truncated 3D projections for four distinct acquisition geometries 
with the following parameters: 

(i) Spherical tomography with sources on a full spherical surface. Ray dis¬ 
cretization: 3 degrees in the polar direction, 5 degrees in the azimuthal 
direction; scanning radius: 400 voxels; number of detector rows: 256; 
source-detector distance = 900 voxels. 

(ii) Spiral tomography with sources on a helix. Helical pitch: 35 voxels; 
8 turns to scan the whole object; number of source positions: 128 per 
complete turn; scanning radius = 384 voxels; number of detector rows 
= 16; source-detector distance: 768 voxels. 

(iii) C-arm tomography with sources on a circle. Scanning radius: 1472 
voxels; number of source positions: 360; detector size: 256 rows, 256 
columns; detector spacing: 1 voxel; source-detector distance: 1472 vox¬ 
els. 
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(iv) Twin circles tomography with sources on two concentric circles, con¬ 
tained in orthogonal planes on M 3 . Common radius of the two circles: 
1472 voxels; number of source positions: 360 per circle; detector size: 
256 rows, 256 columns; detector spacing: 1 voxel; source-detector dis¬ 
tance: 1472 voxels. 

For each one of these four acquisition setups, we selected four concentric 
spherical ROI C with radius values (in voxels) equal to 45, 60, 75, 90. Each 
such ROI C was used to truncate the discretized projection data Y = Df 
to the rays intersecting C and thus to generate a discretized version of the 
ROI-truncated data Yq = Dcf • Note that non truncation corresponded to 
a much larger spherical radius (221 voxels) covering the entire 3D density 
volume. 


6.2 Numerical implementations of our ROI reconstruc¬ 
tion algorithm. 


For each 3D discrete density function /, each cone-beam acquisition setup 
and each spherical ROI C, we have implemented our iterative ROI recon¬ 
struction algorithm to compute a reconstruction Zcf of / using only the 
ROI-truncated data Dcf■ According to our general scheme (21) and Corol¬ 
lary [lj we apply the iterative ROI reconstruction formula f 3+ \ = fo + Ufj 
where U = <jZ(D — Dc ), cr is a regularization operators and Z is the inverse 
of the non-truncated cone-beam transform implemented using the following 
specific methods, according to the acquisition setup. 

For each one of our simulated cone-beam acquisition setups, the inverse 
Z of the non-truncated cone-beam transform was implemented as follows: 


(i) Spherical tomography: Inversion of non-truncated cone-beam transform 
by filtered back-projection (FBP) (3], 

(ii) Spiral tomography: Inversion of non-truncated cone-beam transform by 
Katsevich’s inversion formula [23]. 

(iii) C-arm tomography: Inversion of non-truncated cone-beam transform 
by a well known FDK algorithm [30]. 


(iv) Two circles tomography: Inversion of non-truncated cone-beam trans¬ 
form by a discretized version of Grangeat’s formula outlined in Defrise 
and Clack [31]. 



Choice of a regularization operator a on the euclidean ball B. For 

all our acquisition setups, we used similar discretized versions of the regu¬ 
larization operator a : L 2 (B) —y W 4 (B) based on wavelet thresholding. That 
is, to compute ah for any h in L 2 (B), we first expanded h using standard 
Daubechics wavelets Daub4 [32] in R 3 to generate the wavelet decomposition 
of h 

h = y, a(m, n, i)(f>m,n,i- 

m,n,i 

Each wavelet 4> m , n ,i in this family is a C 4 function indexed by the discretized 
position (m, n) of its compact support and by an integer scale parameter 
i > 0. At the coarsest scale i — 0, no truncation or shrinkage was applied 
to the wavelets coefficients a(m,n,i). At finer scales i > 1, the wavelets 
coefficients a(m,n,i) were set to zero whenever \a(m,n,i)\ < THRj. where 
the thresholds THRi were selected to discard 90% of wavelet coefficients. We 
found that the performance of the algorithm is not very sensitive to the choice 
of the percentage of discarded wavelet coefficients, that is, performance would 
remain essentially the same by discarding 70-90% of wavelet coefficients. 
Additional details can be found in [ 33] . 

The new wavelet expansion generated by this coefficients truncation de¬ 
fined the function ah which obviously belonged to W 4 (B). This operator 
a is non linear but can be well approximated by linearized versions which 
implement smooth shrinkage of the wavelets coefficients instead of abrupt 
truncation (see [33].) 


Stopping rules for our ROI reconstruction Let C be the spherical 
ROI. As a stopping criterion, we adopted a standard rule so that the algo¬ 
rithm (21) is to stop the iteration over the index j when fj and fj + \ become 
close enough within C; in particular, as long as 


ll/j+i _ /jlUhC) < b 

for some small tolerance 6, e.g. b = 0.02. We automatically stop the ROI 
iterative reconstruction at j = 40 to avoid unnecessary computation as we 
found that, for all our numerical experiments, as soon as the radius of C was 
slightly superior to a critical radius, 40 iterations were amply sufficient to 
achieve convergence. 
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6.3 Performances of numerical ROI reconstructions 


For each one of our three densities /, each one of our four x-ray acquisition 
setups, and each one of our selected spherical ROls C, our simulations gener¬ 
ated a discrete version g = D c f of the ROI-truncated cone-beam projections 
of /. Then the numerical application of our iterative ROI reconstruction algo¬ 
rithm to these truncated data g provided a discretized approximation Zcf of 
the ” unknown” /. To assess the accuracy of our discretized ROI reconstruc¬ 
tion Zcf, we have evaluated an ROI Relative L 1 Error of Reconstruction 
within C defined by the following ratio RLi of two discretized T 1 (C) norms 


RLi = 


11/ - Zc/IU.(Q 
ll/lli'to 


For each one of the 48 ROI reconstruction cases indicated above, we have 
recorded the Relative L 1 reconstruction error computed within the ROI in 
Table [I] As indicated above, our iterative algorithm (21) uses different nu¬ 
merical routines to implement the non-truncated inverse operator Z depend¬ 
ing on the acquisitions geometry. Namely, in the case of sources on a sphere, 
Z is implemented using the FBP algorithm; for sources on a spiral curve, 
Z is implemented using the Katsevich’s inversion formula; for sources on 
a circular curve, Z is implemented using the FDK algorithm; for sources 
on a twin-circle curve, Z is implemented using a version of Grangeat’s for¬ 
mula. The number of iterations needed to achieve convergence of our ROI 
reconstruction algorithm was bounded above by 40 but the algorithm was 
found to converge (according to the stopping rule given above) with a much 
smaller number of iterations, typically between 10-12 iterations for sources 
on a curve. 

For each one of the various combinations of density data and acquisition 
setup, these accuracy results yield an estimate of the critical ROI radius p 
enabling a relative ROI reconstruction accuracy inferior or equal to 10%. 
Critical radius estimates are displayed in Table [2] 

The best performances of our ROI reconstruction algorithm naturally 
occur for spherical acquisition. Indeed for the somewhat academic spherical 
setup, the number of projections available is much larger than for the three 
other setups where sources are located on a curve. 

For the twelve situations evaluated here, we obtain a range from 52 to 82 
voxels for the critical radius p yielding a 10% accuracy in ROI reconstruction. 
This compares very favourably to the maximal ROI radius corresponding to 
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Density 

data 

ROI 

radius 

Sources locations 

Spherical 

Spiral 

Circle 

Twin circles 


45 

vox 

10.3% 

10.9% 

13.2% 

14.8% 

Shepp-Logan 

60 

vox 

8.6% 

9.1% 

11.6% 

14.7% 

75 

vox 

7.6% 

8.3% 

7.4% 

8.9% 


90 

vox 

7.3% 

8.0% 

4.4% 

4.8% 


45 

vox 

10.8% 

11.4% 

11.6% 

12.5% 

Mouse tissue 

60 

vox 

8.8% 

9.7% 

11.1% 

9.4% 

75 

vox 

7.9% 

8.8% 

8.4% 

8.3% 


90 

vox 

7.5% 

8.4% 

7.1% 

7.8% 


45 

vox 

11.4% 

11.9% 

12.9% 

15.0% 

Human jaw 

60 

vox 

9.6% 

10.8% 

12.8% 

13.3% 

75 

vox 

9.0% 

9.7% 

10.2% 

10.2% 


90 

vox 

8.2% 

8.5% 

9.8% 

9.8% 


Table 1: Relative L 1 error of ROI reconstruction. The table shows the recon¬ 
struction accuracy within the ROI using four ROI radii for three 3D density 
data and four cone-beam acquisition geometries. Each density data set has 
size 256 3 . 


Density data 

Source locations 

Spherical 

Spiral 

Circle 

Twin circles 

Shepp-Logan 

52 vox 

56 vox 

67 vox 

73 vox 

Mouse tissue 

52 vox 

57 vox 

66 vox 

49 vox 

Human jaw 

57 vox 

70 vox 

82 vox 

82 vox 


Table 2: Critical radius of convergence. For three 3D density data and four 
cone-beam acquisition geometries, the table shows the critical ROI radius 
above which the relative accuracy of the iterative ROI reconstruction was 
found to be less than 0.1. 


31 





non truncation. Indeed when one goes from non truncation to a critical 
spherical ROI, the reduction in irradiated volume ranges from 70% to 98%, 
indicating a quite strong “formal” reduction in x-ray exposure, while the loss 
in relative reconstruction accuracy is only of the order of 7%. 

Note also that the actual critical radius estimates obtained here by simu¬ 
lations are much smaller that the theoretical upper bounds used in the proof 
of Theorem CD 

For a fixed ROI radius, the ROI Relative Reconstruction error is lower for 
the Shepp-Logan phantom than for Mouse Tissue or Human Jaw density data. 
Indeed, when the ROI radius is larger than the critical radius, our iterative 
ROI reconstruction essentially converges within the ROI to a regularization 
erf of /. The ROI reconstruction error in L l (C) can roughly be viewed as the 
sum of two terms, a ‘convergence’ error \\Zcf— c r /||i 1 (c) and a ‘regularization’ 
error \\f — erf\\ L i( C )- To highlight the regularization effect, we have computed 
the relative ‘regularization error’ within C given by 

11/ -o7|| L i( C ) 

II/IUro) 

For an ROI radius of 70 voxels, this regularization error is equal to 1.1% 
for the 3D Shepp-Logan phantom, and to 2.4% for the Mouse Tissue and 
Human Jaw 3D data, because the piecewise constant Shepp-Logan phantom 
density can be approximated by our wavelet-based regularization operator 
much more effectively than the more textured Mouse Tissue and Human Jaw 
densities. So for the Human Jaw data the regularization error contributes 
about half of the relative L 1 reconstruction error. 

To illustrate visually the overall performance of our iterative ROI recon¬ 
struction from truncated cone-beam data, we have include several examples. 
For all these examples, the size of the images is 256 3 voxels and the ROI 
radius is 45 voxels. In Figures [2]{3] we show horizontal, coronal, and sagit¬ 
tal planes from our 3D reconstruction of the Shepp-Logan 3D Phantom and 
Mouse Tissue data using simulated Twin Circle acquisition. We also include 
line profiles to compare our reconstruction against ground truth and one-step 
inversion formula. 

In Figures [4]|5j we show horizontal sections from the reconstructed vol¬ 
umes of the Shepp-Logan 3D Phantom and Mouse Tissue data using simu¬ 
lated spiral and C-arrn acquisitions. 

In all these figures, the comparison of results from our iterative ROI re¬ 
constructions with those obtained by the classical one-step inversion formu- 
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Naive 

reconstruction 


Our approach 


Ground truth 


Line intensity 
profiles 


Figure 2: Visual comparison of ROI reconstruction for 3D Shepp-Logan 
phantom using simulated Twin Circles acquisition and truncation of pro¬ 
jection data. ROI radius = 45 voxels. Middles sections are shown from the 
xy, yz and xz planes. From left to right: inversion by one-step Grangeat for¬ 
mula; our iterative ROI reconstruction; ground truth. The last column shows 
intensity profiles corresponding to the middle row of the images. Green: one- 
step Grangeat formula; blue: our algorithm; red: ground truth. 
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Line intensity 
profiles 


Figure 3: Visual comparison of ROI reconstruction for Mouse Tissue data 
using simulated Twin Circles acquisition and truncation of projection data. 
ROI radius = 45 voxels. Middles sections are shown from the xy , yz and 
xz planes. From left to right: inversion by one-step Grangeat’s formula; our 
iterative ROI reconstruction; ground truth. The last column shows intensity 
profiles corresponding to the middle row of the images. Green: one-step 
Grangeat formula; blue: our algorithm; red: ground truth. 
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las originally devised for reconstruction from non-truncated cone-beam data 
show that, as expected, the one-step inversion formulas for non-truncated 
data perform poorly when applied to ROI-truncated cone-beam data, and 
display multiple visual artifacts especially near the ROl boundary. By con¬ 
trast, our iterative ROI reconstruction results are very satisfactory even for 
relatively small ROI radii. Compared to the ground truth, our ROI recon¬ 
struction shows some blurring which is due to the wavelet-based regulariza¬ 
tion step. Since our wavelet filters have finite support and length 4 pixels, 
they have a rather limited impact on spatial resolution. The blurring effect 
is consistent with our theoretical prediction since our algorithm generates an 
approximation of the exact solution which is a smoother version of the true 
image. 

Note also that the images reconstructed using one-step Katsevich formu¬ 
las in Figure [4] exhibit streak artifacts. This is a common and known issue 
in spiral tomography, cf. [34]. Our regularized reconstruction significantly 
reduces these artifacts through the wavelet-based regularization. 

7 Conclusion 

In this paper, we have examined the problem of ROI tomographic recon¬ 
struction using truncated cone-beam data, a problem of high relevance in 
many applications. For both our theoretical and numerical analysis, we con¬ 
sidered fairly generic cone-beam acquisition setups, with sources located on 
arbitrary bounded smooth curves T in M 3 verifying classical Tuy’s condition. 
In all these cases, it is known that the non-trucated cone-beam transform 
Df of smooth densities / admits an explicit inverse Z but Z cannot directly 
reconstruct / from ROI-truncated data. 

To deal with the reconstruction from ROI-truncated data, we have de¬ 
veloped and rigorously analyzed a new iterative ROI reconstruction method 
valid for densities / in where B is a bounded ball in M 3 , which it¬ 

erates a linear contraction endomorphism U of L°°(B). The operator U is 
constructed by combining forward ROI-truncated projections, backward in¬ 
version by the operator Z and appropriate regularization operators defined 
in image and/or projection space. Our main theoretical result is that, given 
e > 0, for spherical regions of interest C C B with radius larger than a criti¬ 
cal radius p(e): (i) our iterative ROI reconstruction from ROI-truncated data 
converges in L°°(B ) to a density estimate / such that \\f — f\\oo < e ||/||oo; 
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Shepp-Logan 

phantom 


Mouse tissue 


Naive 

reconstruction 


Our approach Ground truth 



Figure 4: Visual comparison of ROI reconstruction for 3D Shepp-Logan 
phantom and mouse tissue using simulated spiral acquisition and trunca¬ 
tion of projection data. A representative horizontal section from the 3D 
reconstructed volume is shown. From left to right: inversion by one-step 
Katsevich formula; our iterative ROI reconstruction; ground truth. 
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Naive 

reconstruction 


Our approach Ground truth 


Figure 5: Visual comparison of ROI reconstruction for 3D Shepp-Logan 
phantom and mouse tissue using simulated C-arm acquisition and trunca¬ 
tion of projection data. A representative horizontal section from the 3D 
reconstructed volume is shown. From left to right: inversion by one-step 
FDK algorithm; our iterative ROI reconstruction; ground truth. 
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(ii) our iterative ROI reconstruction algorithm generates a bounded linear 
operator Zc '■ L°°(1Zb) —>■ L°°(B), where IZb is the Riemannian manifold of 
all x-rays emitted by sources on a curve T outside B. The operator Zc is an e- 
inverse of the ROI-truncated cone-beam transform Dq : L°°(B) — > L°°(1Zb ), 
that is ||J — Z c D c \\l°°(b) < £■ These results also extend to the case of 
spherical acquisition in M 3 and to the ray transform. 

Even though iterative methods for ROI CT reconstruction already ap¬ 
peared in the literature, up to the knowledge of the authors, no theoretical 
result was known so far about the existence of a critical radius ensuring the 
convergence of an iterative ROI CT reconstruction scheme. 

We numerically verified our theoretical results using simulated 3D ac¬ 
quisition of ROI-truncated cone-beam data for four classical acquisition ge¬ 
ometries (spherical, spiral, circular arm, twin orthogonal circles), using three 
different density functions and multiple ROI radii and locations. All numer¬ 
ical experiments show that, for e moderately small, e.g., e = 0.1, the critical 
ROI radius p(e) is fairly small with respect to the support of the density 
function. 
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Appendix: Sobolev imbeddings 

In Section |3.1[ to show the regularity of the linear operator Z, we make use 
of Sobolev imbedding theorems. We quote a special case of a result by Aubin 
(Theorem 2.34 in [35J) in this context. 

Theorem 2 ([35]). If A4 is a compact Riemannian manifold of dimension 
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n with C 1 -boundary and interior M., then 


W k (M ) C C a (M) 

and this imbedding is compact if k — a > n/2. 

By the compactness of the embedding, it is also continuous, meaning that 
if k — a > n/2 then there exists c > 0 such that for each / e W A (M.), 

\\f\\c a (M) ^ c \\f\\w 4 (M) ■ 

In particular, if n — 3 and k = 4, then we can choose a = 2 and have the 
following result. 

Corollary 2. If Ai is a compact Riemannian manifold of dimension n with 
C 1 -boundary and interior M. then 

\\f\\c 2 (M) — c \\f\\w 4 {M ) • 
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